1 Local MDS for Zero Digits

As in Assignment 2, we will consider only the zeros in the ZIP data set.

# Utility function to plot images. 
plot.zip <- function(x,use.first=FALSE, ...){
  x<-as.numeric(x)
  if (use.first){
    x.mat <- matrix(x,16,16)
  }else{
    x.mat <- matrix(x[-1],16,16)
  }
  image(1:16,1:16,x.mat[,16:1],
        col=gray(seq(1,0,l=12)))
  if (!use.first) title(x[1])
}

zip.train <- read.table("../../unsup/zip.train")
zip.zeros <- zip.train %>% dplyr::filter(V1 == 0) %>% dplyr::select(-1)
dim(zip.zeros)
#> [1] 1194  256
delta <- dist(zip.zeros)

1.1 2-dimensional Configuration of the Data

Using parameters \(k = 5\) and \(\tau = 0.05\) in the function stops::lmds, a 2-dimensional configuration of the data is found.

q <- 2
k <- 5
tau <- 0.05
lmbds.res <- lmds(as.matrix(delta), k = k, tau = tau, ndim = q, itmax = 700, verbose = 10)
#> [1] "niter=100 stress=-251754.80998"
#> [1] "niter=200 stress=-265539.92165"
#> [1] "niter=300 stress=-269891.73297"
#> [1] "niter=400 stress=-271730.14941"
#> [1] "niter=500 stress=-272742.21209"
#> [1] "niter=600 stress=-273206.77601"
#> [1] "niter=700 stress=-273634.96777"
conf <- lmbds.res$conf

The scatterplot is shown below. 9 points are selected in order to try to “cover” the variability in the scatter plot. The points we have chosen are shown in red. They are chosen from left to right along the first dimension and from bottom to top along the second dimension.

# Code to select points. 
find.nearest <- function(new.pos, data){
  # Compute distance to points and select nearest index.
  return(which.min(colSums((t(data) - new.pos)^2)))
}


chosen.indices.localMDS <- c(find.nearest(c(-100,0), conf), 
                                find.nearest(c(-40,0), conf), 
                                find.nearest(c(0,0), conf), 
                                find.nearest(c(20,0), conf), 
                                find.nearest(c(60,0), conf), 
                                find.nearest(c(40,-25), conf), 
                                find.nearest(c(40,0), conf), 
                                find.nearest(c(40,20), conf),
                                find.nearest(c(40,50), conf))


chosen.points.localMDS <- data.matrix(conf[chosen.indices.localMDS,])

plot(conf, as = 1, main=paste0("Local MDS, k=",k,", tau=",tau), xlab = "Dim1", ylab = "Dim2")
points(chosen.points.localMDS, col = 2, pch = 16)

The ZERO digits corresponding to each of the points are plotted below.

par(mfrow = c(3,3))
apply(zip.zeros[chosen.indices.localMDS,],1,plot.zip, use.first = T)

#> NULL

The first five zeros (row-wise) travel from left to right in the first dimension of the scatter plot gained from local MDS. It seems like the first coordinate describes zeros from narrow to wide (in the west-east direction of the drawings). The last four zeros that are plotted show the movement from small values to larger values in the second dimension. It seems like the second coordinates describes zeros from thin thick (in the south-north direction of the drawings). Thus, the first dimensions describe the changing sizes of the zeros, while the second dimensions describes the different line thicknesses of the zeros. This can be seen more clearly when plotting more points, but we have only shown 9 here.

1.1.1 (OPTIONAL) Relate the results from Local MDS with those obtained by the first 3 principal components. In particular, could you represent in any way the results obtained by Local MDS in the 3-dimensional scatterplot of (PC1, PC2, PC3)?

Below the dimensions found using Local MDS are plotted against the first three principal components. We observe that the first dimension of the Local MDS posetively correlates with the first PC. Similarly the second dimension of the Local MDS correlates negatively with the second PC. Third PC does not show as clear of a correlations, but it absolutely show a pattern. To create the plots we fit a GAM using splines for each principal component against the coordinates in the Local MDS dimensions. Finally we plot all three together in the PC1,PC2,PC3 plane. The plot seems reasonable, but is hard to read and interpret.

zip.PC <- princomp(zip.zeros[,-1])
scores <- zip.PC$scores[,1:3]

smooth.pc1 <- gam(scores[,1] ~ s(conf[,1],conf[,2]))
smooth.pc2 <- gam(scores[,2] ~ s(conf[,1],conf[,2]))
smooth.pc3 <- gam(scores[,3] ~ s(conf[,1],conf[,2]))

s.pc1 <- smooth.pc1$fitted.values
s.pc2 <- smooth.pc2$fitted.values
s.pc3 <- smooth.pc3$fitted.values

points3D(conf[,1], conf[,2], s.pc1,
         xlab = "Dim1", ylab = "Dim2", zlab ="PC1")

points3D(conf[,1], conf[,2], s.pc2,
         xlab = "Dim1", ylab = "Dim2", zlab ="PC2")

points3D(conf[,1], conf[,2], s.pc3,
         xlab = "Dim1", ylab = "Dim2", zlab ="PC3")

points3D(scores[,1], scores[,2], scores[,3], pch=19,cex=.4, col = 1, 
       xlim = range(scores[,1]), ylim = range(scores[,2]), zlim = range(scores[,3]),
       xlab = "PC1", ylab = "PC2", zlab ="PC3")
points3D(s.pc1, s.pc2, s.pc3, col = 2, add = T, alpha = 0.5)

2 ISOMAP for Zero Digits

The same tasks are solved with ISOMAP instead, in order to compare to the results obtained when using local MDS.

k.iso <- 5
ismp <- isomap(delta, k = k.iso, ndim = q)
ismp.points <- ismp$points

The two-dimensional configuration found using ISOMAP can be seen below.

aux.plot <- plot(ismp,n.col=3,main="Output of ISOMAP Algorithm")
points(aux.plot,"sites",pch=19,cex=.6)

chosen.indices.ISO <- c(find.nearest(c(-45,20), ismp.points), 
                                find.nearest(c(0,20), ismp.points), 
                                find.nearest(c(45,20), ismp.points), 
                                find.nearest(c(-45,-10), ismp.points), 
                                find.nearest(c(0,-10), ismp.points), 
                                find.nearest(c(45,-10), ismp.points), 
                                find.nearest(c(-45,-45), ismp.points), 
                                find.nearest(c(0,-45), ismp.points),
                                find.nearest(c(45,-45), ismp.points))
chosen.indices.ISO <- c(find.nearest(c(-100,0), ismp.points), 
                                find.nearest(c(-40,0), ismp.points), 
                                find.nearest(c(-10,0), ismp.points), 
                                find.nearest(c(10,0), ismp.points), 
                                find.nearest(c(40,0), ismp.points), 
                                find.nearest(c(25,-45), ismp.points), 
                                find.nearest(c(25,-20), ismp.points), 
                                find.nearest(c(25,0), ismp.points),
                                find.nearest(c(25,25), ismp.points))

chosen.points.ISO<- data.matrix(ismp.points[chosen.indices.ISO,])

plot(ismp.points, main=paste0("ISOMAP, k=",k.iso), xlab = "Dim1", ylab = "Dim2")
points(chosen.points.ISO, col = 2, pch = 16)

After selecting points the red points shown above, the corresponding zeros are plotted in a similar fashion as for local MDS. Notice that the red points have now been selected in a grid-like way. Starting from the top left followed by the top middle, top right and then middle left and so on. Since the points are less spread to the left this means that the three left selected points are more close to each other.

par(mfrow = c(3,3))
apply(zip.zeros[chosen.indices.ISO, ],1,plot.zip, use.first = T)

#> NULL
par(mfrow = c(1,1))

#As earlier, the first five images show the traversion from west to east in the first dimension, while the last four images show the traversion from south to north in the second dimension. The first dimension seems to describe the … and the second dimension seems to describe the …

As before, the zeros corresponding to the red points are plotted. The grid of zeros now directly corresponds to the, somewhat skewed, grid of red points above. When inspecting, we can see that the traversion from west to east seems to describe the width of the zeros while the traversion from south to north are probably describing the hight of the zero.

TOLKNINGEN HER MANGLER!

3 Select the Tuning Parameters for Zero Digits

3.1 ISOMAP

The local continuity meta criteria is used to select the tuning parameter \(k\) in ISOMAP. The first cell of code is used to calculate the adjusted local continuity meta criteria (function taken from the lecture material).

LCMC <- function(D1,D2,Kp){
  D1 <- as.matrix(D1)
  D2 <- as.matrix(D2)
  n <- dim(D1)[1]
  N.Kp.i <- numeric(n)
  for (i in 1:n){
    N1.i <- sort.int(D1[i,],index.return = TRUE)$ix[1:Kp]
    N2.i <- sort.int(D2[i,],index.return = TRUE)$ix[1:Kp]
    N.Kp.i[i] <- length(intersect(N1.i, N2.i))
  }
  N.Kp<-mean(N.Kp.i)
  M.Kp.adj <- N.Kp/Kp - Kp/(n-1)
  
  return(list(N.Kp.i=N.Kp.i, M.Kp.adj=M.Kp.adj))
}

Next, the function is used to find optimal \(k\). We fix the dimension \(q = 2\). The search for the optimal \(k\) is restricted to integers between 3 and 10. Note that we begin at \(k = 3\) because a lower \(k\) gives fragmented data in the ISOMAP.

Kp <- 5
k.search <- 3:10
nk <- length(k.search)
LC <- rep(NA, nk)
ISOMAP.k <- vector("list",nk)

for (i in 1:nk){
  ISOMAP.k[[i]] <- isomap(delta, ndim=q, 
                            k = k.search[i])
  D2.k <- dist(ISOMAP.k[[i]]$points[,1:q])
  LC[i] <- LCMC(delta,D2.k,Kp)$M.Kp.adj
}

i.max <- which.max(LC)
k.max <- k.search[i.max]
ISOMAP.max <- ISOMAP.k[[i.max]]

plot(k.search, LC, type="b", main=paste0("k.max=",round(k.max,4)))
abline(v=k.max,col=2)

The plot above shows that the optimal \(k\) is \(k =\) 3. The graphical description of the low dimensional configuration corresponding to the optimal \(k\) is given in the following.

ismp.opt <- isomap(delta, k = k.max, ndim = q)
ismp.opt.points <- ismp.opt$points

aux.plot.opt <- plot(ismp.opt,n.col=3,main=paste0("Output of ISOMAP Algorithm with Optimal k=",k.max))
points(aux.plot.opt,"sites",pch=19,cex=.6)

A similar description for 9 chosen digits can be done here. HER MANGLER DET NOE

3.2 Local MDS

The local continuity meta criteria can also be used to select the optimal parameters \(\tau\) and \(k\) in Local MDS. We restrict the search of \(k\) to the integers \(5, 10, 15\) and the search for \(\tau\) is restricted to \(0.1, 0.5, 1\).

Kp <- 10
K <- c(5,10,15)
tau <- c(.1,.5,1)

LC <- matrix(0,nrow=length(K),ncol=length(tau))
LocalMDS.k.tau <- array(vector("list",1),dim=dim(LC))

for (i in 1:length(K)){
  for (j in 1:length(tau)){
    LocalMDS.k.tau[[i,j]] <- lmds(as.matrix(delta), k = i, tau = j, ndim = q, itmax = 100)$conf
    D2.k.tau <- dist(LocalMDS.k.tau[[i,j]])
    LC[i,j] <- LCMC(delta,D2.k.tau,Kp)$M.Kp.adj
  }
}

ij.max <- arrayInd(which.max(LC),.dim=dim(LC))
k.max <- K[ij.max[1]] 
tau.max <- tau[ij.max[2]] 
LocalMDS.max <- LocalMDS.k.tau[[ij.max[1],ij.max[2]]]

print(paste0("k.max=",k.max,"; tau.max=",tau.max))
#> [1] "k.max=5; tau.max=0.5"

HER MANGLER DET NOE

4 t-SNE for Zero Digits.

The same tasks are solved with t-SNE, in order to compare to the results obtained earlier.

tsne <- Rtsne(zip.zeros, dims = q, perplexity = 40, theta = 0)
tsne.points <- tsne$Y
plot(tsne.points, main="t-SNE", xlab = "Dim1", ylab = "Dim2")

plot(tsne.points, as = 1, main="t-SNE", xlab = "Dim1", ylab = "Dim2")

We choose some points in order to “cover” the variance in the two-dimensional representation computed by t-SNE.

chosen.indices.tsne <- c(find.nearest(c(-60,0), tsne.points), 
                                find.nearest(c(-30,0), tsne.points), 
                                find.nearest(c(-10,0), tsne.points), 
                                find.nearest(c(20,0), tsne.points), 
                                find.nearest(c(60,0), tsne.points), 
                                find.nearest(c(-30,-30), tsne.points), 
                                find.nearest(c(-30,0), tsne.points), 
                                find.nearest(c(-30,15), tsne.points),
                                find.nearest(c(-30,40), tsne.points))

chosen.points.tsne <- data.matrix(tsne.points[chosen.indices.tsne,])

plot(tsne.points, main="t-SNE", xlab = "Dim1", ylab = "Dim2")
points(chosen.points.tsne, col = 2, pch = 16)

The corresponding zeros to the selected points are plotted below.

par(mfrow = c(3,3))
apply(zip.zeros[chosen.indices.tsne, ],1,plot.zip, use.first = T)

#> NULL
par(mfrow = c(1,3))

As earlier, the first five images show the traversion from west to east in the first dimension, while the last four images show the traversion from south to north in the second dimension. The first dimension seems to describe the … and the second dimension seems to describe the …

5 Select Tuning Parameter for Zero Digits

The local continuity meta criteria is used to select the tuning parameter perplexity in t-SNE. We are still using the fixed lower dimension \(q = 2\) and we are setting theta \(= 0\). The search for the optimal value of perplexity is constrained to the values c(20, 40, 60).

perplex <- c(20,40,60)
Kp <- 10

LC <- numeric(length(perplex))
Rtsne.k <- vector("list",length(perplex))

for (i in 1:length(perplex)){
    Rtsne.k[[i]] <- Rtsne(zip.zeros, perplexity=perplex[i], dims=q,
                          theta=0, pca=FALSE, max_iter = 500)
    D2.k <- dist(Rtsne.k[[i]]$Y)
    LC[i] <- LCMC(delta,D2.k,Kp)$M.Kp.adj
}

i.max <- which.max(LC)
perplexity.max <- perplex[i.max[1]] 
Rtsne.max <- Rtsne.k[[i.max]]

plot(perplex,LC, main=paste0("perplexity.max=",perplexity.max))
abline(v=perplex[i.max],col=2)

HER MANGLER DET NOE